Optimal waveform estimation for classical and quantum systems via time-symmetric 
smoothing. II. Applications to atomic magnetometry and Hardy's paradox 
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The quantu m smo othing theory [Tsang, Phys. Rev. Lett. 102, 250403 (2009); Phys. Rev. A, in 
press (e-print arXiv:0906.4133 )] is extended to account for discrete jumps in the classical random 
process to be estimated, discrete variables in the quantum system, such as spin, angular momentum, 
and photon number, and Poissonian measurements, such as photon counting. The extended theory 
is used to model atomic magnetometers and study Hardy's paradox in phase space. In the phase- 
space picture of Hardy's proposed experiment, the negativity of the predictive Wigner distribution 
is identified as the culprit of the disagreement between classical reasoning and quantum mechanics. 
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I. INTRODUCTION 

In previous papers [J, Q j I have proposed a quantum 
smoothing theory, which can be used to optimally esti- 
mate classical signals coupled to quantum sensors under 
continuous measurements, such as gravitational wave de- 
tectors and atomic magnetometers. Smoothing can be 
significantly more accurate than current quantum filter- 
ing methods [1, i, i, d, 0, S, S 0, [HI when the classical 
signal is a stochastic process and delay is permitted in the 
estimation. While Refs. P, 0] focus on diffusive classical 
random processes, quantum systems with continuous de- 
grees of freedom, and Gaussian measurements, the aim of 
this paper is to extend the theory to account for discrete 
variables in the systems and the measurements. In par- 
ticular, I shall consider discrete jumps in the classical ran- 
dom process, discrete variables in the quantum system, 
such as spin, angular momentum, and photon number, 
and Poissonian measurements, such as photon counting. 
Such extensions are especially important for the model- 
ing of atomic magnetometry [1' B E [IH, [H, [11 [11 ■ 

In the case of atomic magnetometry, the importance of 
estimation delay was discovered by Petersen and M0lmer 
[3] , who found that the estimation of a fluctuating mag- 
netic field modeled as an Ornstein-Uhlenbeck process be- 
comes more accurate when the estimation is delayed and 
observations at later times are taken into account. I shall 
generalize their results using the quantum smoothing the- 
ory, derive the optimal strategy of delayed estimation for 
atomic magnetometry, and discuss practical methods of 
implementing the strategy. 

A different kind of estimation problem comes up in 
Hardy's paradox [l5| . in which one estimates the posi- 
tions of an electron and a positron in interferometers 
based on posterior measurement outcomes and obtains 
paradoxical results. I shall demonstrate that the salient 
features of the paradox can be reproduced mathemat- 



ically using the quantum smoothing theory in discrete 
phase space, which is arguably the most natural way 
of modeling classical properties of quantum objects. It 
is shown that the negativity of the predictive Wigner 
distribution can be regarded as the culprit of the dis- 
agreement between classical reasoning and quantum me- 
chanics. This phase-space approach is somewhat differ- 
ent from Aharonov et aZ.'s weak value approach [l6j . 
Whether the two can be reconciled remains to be seen. 

This paper is organized as follows: Section HIl reviews 
the classical filtering and smoothing equations when the 
system process has jumps and the observations have Pois- 
sonian statistics, as derived by Snyder [13, [ll] and Par- 
doux . Section IIIII generalizes such equations to the 
quantum regime for smoothing of classical random pro- 
cesses coupled to quantum systems. Sec. IIVI converts the 
quantum equations to equivalent phase-space equations 
for discrete Wigner distributions. Sec. |V] studies the ap- 
plication of the theory to atomic magnetometry. Sec. IVII 
studies Hardy's paradox using quantum smoothing in dis- 
crete phase space. 



II. CLASSICAL FILTERING AND SMOOTHING 
FOR POISSONIAN OBSERVATIONS 

Define Xt as the classical system random process, the 
a priori probability density of which satisfies the differ- 
ential Chapman-Kolmogorov equation [20j 
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where J{x\x' ,t) is the probabihty density per unit time 
that Xt will jump from x' to x. For an obseration process 
with Poissonian noise, the conditional probability density 
is 



P{5N^,t\xt) = exp [-\^{xt,t)5t] 



with the continuous-time limit 



[\^{xt,t)5t\'^- 



, (2.3) 



dN% = dN^u (2.4) 
P{dN^t = 0|xt) = 1 - \{xut)dt, (2.5) 
P{dN^t = l\xt) = X^ixt,t)dt. (2.6) 

Defining the observation record in the time interval < 
t <t as 



dN[t^^t^={dNutQ<t<t} 



(2.7) 



and the filtering probability density as the probability 
density of xt conditioned upon past observations, given 
by 



F{x,t) = P{xt^x\dN[t,,t)), 



(2.8) 



the Ito stochastic differential equation for F{x, t) is called 
the Snyder equation and given by [13, fisj 

dF = dtCcF + J2 {dN^t - dt {X^)p) (A^);' 

X (A^-(A,,)^)F, (2.9) 



= / dxX^,{x,t)F{x,t). 



(2.10) 



A linear equation for an unnormalized F(x, t) was derived 
by Pardoux and given by 

dj = dtCcf + ^ {dN^t - dt) (A^ - 1)/, (2.11) 



with 



F{x,t) 



J dxf{x,t)' 



(2.12) 



To perform smoothing in the time-symmetric form fl9'|, 
first solve for an unnormalized retrodictive likelihood 
function P{dN[f rp'^\xt — x) on g{x,t) using the adjoint 
of Eq. (EH]), 

~dg = dtC*cg + {dN^t - dt) (A^ - l).g, (2.13) 



to be solved backward in time with final condition 
g{x, T) oc 1. The smoothing probability density at time r 
given the observation record dN^j-^ rp^ in the time interval 
to <T <T is then 

vf i^AT \ 9{x,T)f{x,T) 

P[Xt = x\dN[to.T)) r j^^/^ (2-14) 



J dxg{x,T)f{x,T) ' 



III. HYBRID CLASSICAL-QUANTUM 
FILTERING AND SMOOTHING FOR 
POISSONIAN OBSERVATIONS 

Using the same approach as Refs. [H, 01, it is not dif- 
ficult to generalize the above classical equations to the 
quantum regime for hybrid classical-quantum filtering 
and smoothing. Define Xt as the classical system process 
that one wishes to estimate, which is coupled to a quan- 
tum system under measurements. As before, the quan- 
tum backaction from the quantum system to the classical 
one is assumed to be negligible. Define the hybrid density 
operator that describes the joint statistics of the classi- 
cal and quantum systems [2l| as p{x,t). The a priori 
evolution of p{xt,t) is governed by 



dp{x,t) 

dt 
Cp{x,t) 



Cpix,t), (3.1) 
Cop{x, t) + Ci{x)p{x, t) + Ccp{x, t), (3.2) 



where Cq is the superoperator that governs the evolution 
of the quantum system, Ci is the superoperator that de- 
scribes the coupling of the classical system to the quan- 
tum system, via an interaction Hamiltonian for example, 
and Cq is the Chapman-Kolmogorov operator defined by 
Eq. (|2.2p . The measurement, on the other hand, is de- 
scribed by the quantum Bayes theorem, 



p{xt\6N^t) 



MiSN^t\xt)pixt)MHSN^t\xt) 
J dxt tr (numerator) ' 



(3.3) 



where the measurement operator with Poissonian statis- 
tics is 



Af((5A^^t|xt) = exp 



-7:I^i{xt,t)Lf,{xt,t)St 



[L^ixt,t)VSt]'^-^ 



(3.4) 



where L^{xt,t) is a hybrid operator, an annihilation op- 
erator for example, and can also depend on Xf. In the 
continuous- time limit, 

M{dNf,t = 0\xt) = i - lLUxt,t)Lf,{xt,t)dt, (3.5) 



(3.6) 



M{dN^,t = l\xt) = L^{xt,t)Vdt. 



After some algebra, the stochastic differential equation 
for the filtering hybrid density operator, defined as 



Fix,t) 
is given by [1, [13] 

dF ^ dtCF + dtY^ 



p{xt = x\dN^to.t)), 



(3.7) 



E 



dN^t - dt 



\lIl,f^\flIl, 



(l,flI-{lIl,)^f), 



(3.8) 
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where 



— I dxii 



LUx,t)L^{x,t)F{x,t) . (3.9) 



Equation (|3.8p is a quantum generalization of the Sny- 
der equation [Eq. (|2.9p ]. A hnear version of Eq. (|3.8p . 
analogous to Eq. (|2.11[) . may be written as [l^j 



df = dtCf + dtY, (lJlI - h^Lj - \fLlL, 



2 M 2 
^ {dN^t - dt) (lJlI - /) , (3.10) 



Fix,t) = 



J dxtY[f{x,t)] 



(3.11) 



The classical incoherent limit of Eq. (I3.10|) is obviously 
Eq. and Eq. (pHT)) can be verified against Eq. 

by normalizing the former using Ito rule. 

To perform smoothing, one also needs to solve for the 
unnormalizcd hybrid effect operator E(dN\f^j'^\xt = x) (x 
g{x,t) using the adjoint of Eq. (PrU]) [l|,[|],' 



dg = dtC*g + dt^ [LlgL^ - ^qLIl,, - ^lIl,^ 



+ Y,idN,t-dt) (LlgL^-g) 



(3.12) 



where the final condition is g{x^ T) oc 1 and the adjoint is 
defined with respect to the Hilbert-Schmidt inner prod- 
uct 



{g{x,t),j{x,t)) = j dxtr [g{x,t)f{x,t)\ , (3.13) 
5(x, t), t)) = (C*g{x, t)Jix, t)) . (3.14) 



The smoothing probability density is then 



filtering equation for F(x, t) is 

dF = dtCF + dtJ2 (l^^FLI - UlL^F - \flIL, 



? (2C'^R-^FC^ - C^^R-'CF - FC^^R'^C 



J2(dN,,~dt{LlL,) ) (LlL 



X ^L,FLl-[LlL,).F 



{C-(C)p^ R^^drjtF + U. 



dt 



dr^t=dyt~-(^C + C^)^, 



(3.16) 
(3.17) 



where C — C{x,t) is a vector of hybrid operators, 
R — R{t) is a positive-definite matrix, dr]t is a vec- 
toral Wiener increment with covariance matrix R(t)dt, 
and H. c. denotes Hermitian conjugate. 
The equation for f{x,t) is 

df = dtCf + dtJ2 [lJlI - ^lIlJ - \fLlL, 

+ y {2C'^R-^fC^ - C''^R-^Cf - fC^'^R-^C^ 
+ J2{dN^,-dt) {LjLl-f) 



+ - (^C^ R-'dytf + n.c. 
and for g{x, t), 

-dg = dt£*g + dtY, UlgL^ - igij^L 
dt 



(3.18) 



2^11 ^,g 



I (2C't^i?-ig(7 - gC^^R~^C - C^^R-^Cg) 
{dN^t - dt) (LlgL^ - g) 



-{gC R-'dyt + U.c. 



(3.19) 



IV. QUANTUM SMOOTHING IN PHASE 
SPACE 



Hx.t) = P{Xr = x\dN[to,T)) 



tr[.g(a:,T)/(a:,T)] 

J dx ti[g{x,T)f{x,T)] 

(3.15) 



Incorporating the Gaussian measurements considered in 
Refs. P, Q] into the equations above is straightforward. 
This is useful, for example, when both photon count- 
ing and homodyne detection are performed in a quan- 
tum optics experiment ^|. With Poissonian observa- 
tions dNt and Gaussian observations dyt, the resulting 



One method of solving Eqs. ([315]) . (|3T8ll . and (|3T9l) 
for hybrid smoothing is to use Wigner distributions ^, . 
For a quantum system with a discrete degree of freedom, 
such as spin, angular momentum, or an A'^-level system, 
one may define the discrete Wigner distribution, accord- 
ing to Feynman and Wootters [1^, as 



f{q,p,x,t) = ^tr f{x,t)W{q,p) 

ge {0,l,...,iV-l}, 
pe{0,l,...,7V-l}. 



(4.1) 

(4.2) 
(4.3) 
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The operator W(q,p) for prime N is 
W{q,p) 

\ + (-l)Pa, + (-l)«+f<7y + i] , TV - 2; 

Egi,g2'52g,gi+g.exp[^p(gi-g2)] |gi>(92|, > 2. 

(4.4) 

(T^;, (Ty, and (Tz are Pauli matrices, and \q2) are eigen- 
states of q, and modular arithmetic with modulus iV is 
implicitly assumed. For a nonprime N ^ the system can 
be decomposed into subsystems with prime iV's and the 
Wigner distribution can be defined using W{q,-p) for each 
subsystem psj . 

An alternative definition in a 2N x 2iV phase space, 
first suggested by Hannay and Berry is 



f{q,p,x,t) 



2N 



■tr 



f{x,t)w{q,p) 



1 



,iV- - 
2 



N 1 iV 



1,. 



N 



(4.5) 



where the matrix elements with noninteger indices are 
assumed to be zero. One may also use either Wigner 
function to describe the energy level n and phase of a 
harmonic oscillator by letting n — q, (/) = 2t:p/N, and 
taking the N ^ oo limit at the end of a calculation [g^, 

m. 

Both definitions have a particularly desirable property 
for the purpose of smoothing, namely. 



tr 



g{x,t)f{x,t) 



N^g{q,p,x,t)f{q,p,x,t), (4.6) 
2iV^.g(<Z,p,x,t)/(<7,p,x,0, (4.7) 



9,P 



SO the smoothing probability density can be written in 
terms of the Wigner distributions as 



or 



h{x, t) 



h{x, t) 



T,q^p9{Q,P,x,T)f{q,p,x,T) 



I J2q.,p 9iQ, X, T)f{q,p, X, t) 

J2q,p 9{<liP, a;, T)f{q,p, X, t) 
I J2q,p 9{q, P, X, T)f{q, p, X, t) 



(4.8) 



(4.9) 



Equations (|4.8p and (|4.9p become equivalent to the clas- 
sical smoothing density given by Eq. (|2.14p . with the 
quantum degrees of freedom marginalized, if / and g or 
/ and g are nonnegative and can be regarded as classi- 
cal probability densities. The hybrid smoothing problem 



can then be solved using classical filtering and smoothing 
techniques. 

If one desires to obtain smoothing estimates of the 
quantum degrees of freedom, a quantum smoothing 
quasiprobability distribution may be defined as 



Hq,P,x,T) 



or 



h{q,p,x,T) 



g{q,p,x,T)f{q,p,x,T) 
I dx Y.q,p 9{<1,P, X, T)f{q, p, X, t) 



g{q,p,x,T)J{q,p,x,T) 
I dx J2q,p gil^P, X, T)f{q,p, X, t) 



(4.10) 



(4.11) 



From the perspective of estimation theory, these defini- 
tions of quantum smoothing distributions are arguably 
the most natural, for they both give the correct classi- 
cal smoothing distribution when the quantum degrees of 
freedom are marginalized, are equivalent to the smooth- 
ing distributions in classical smoothing theory when / 
and g or f and g are nonnegative, and are explicitly nor- 
malized. 

There are many other equally qualified definitions of 
the Wigner distribution in discrete or periodic phase 
space [23, H^l ■ Choosing which definition to use depends 
on the application. The Feynman-Wootters distribution 
is defined only on the eigenvalues of q and p, so it appears 
more physical, but the Hannay-Berry definition is easier 
to calculate analytically for arbitrary N and, as shown in 
Appendix [XI naturally arises from the statistics of weak 
measurements. 



V. ATOMIC MAGNETOMETRY 

A. Optimal smoothing 

An important application of quantum estimation the- 
ory is atomic magnetometry [l, [T H . [ill [l^ [isl. [l^ . 
Consider the setup described in Refs. [3, [lO, [U, [I4I and 
depicted in Fig. [T] An atomic spin ensemble is initially 
prepared in a coherent state with the mean collective 
spin vector along the x axis. Let the magnetic field be 
polarized along y axis and given by 



bt = xit, 



(5.1) 



a component of the classical system process to be esti- 
mated. The magnetic field introduces Larmor precession 
to the spin via the interaction Hamiltonian 



Hi{x) = 
Ci{x)F{x,t) = 



-jbSy 

i 



(5.2) 



Hi{x),F{x,t) 



'-lb 



Sy,F 



(5.3) 



where Sy is the y component of the spin vector operator 
and 7 is the gyromagnetic ratio. Under continuous op- 
tical polarimetry measurements, the stochastic equation 
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Laser Beam 

FIG. 1: (color online). Left: basic setup of atomic magne- 
tometer as considered in Refs. 0, [lol . [ill . H^ . Right: the 
spherical phase space for spin. 



for the filtering density operator F{x, t) has been derived 
by Bouten et al. and is given by 



dF ^ dt^CcF + j-b Sy,F 

|^cos(KTO)i^cos(Km) + sin(Km)i^sin(Km) — F 
^ (dN,,-dt{LlL 



} 



which is in the form of Eq. (|3.8p . with 

m = — ^, L± = —j= [cos(kto) ± sin(KTO)] 
V 2 



(5.4) 



(5.5) 



K is the light-spin coupling parameter and a is the nor- 
malized optical envelope. The linear predictive and retro- 
dictive equations for f{x,t) and g{x,t) become 

df = dt[Ccf+'-^b[Syj 

X l^cos(Krh)/ cos(kto) + sin(KTO)/ sin(Km) — / | 
+ {dN^t-dt)(LjL\-f 



(5.6) 



-dg — dt\ji*Q 



«7, 



Sy,9 



dtl 



X COS 



(Km)^cos(Km) + sin(K7fi)g sin(Km) — g] | 



+ Yl {dN^t~dt)(LlgL^-cj) 



(5.7) 



After solving Eq. (|5.6p forward in time for f{x,T) and 
Eq. (|5.7p backward in time for g{x,T), the smoothing 
probabihty distribution is given by 



h{x, t) = 



ta-[g(a^,T)/(x,T)] 
J dxti[g{x,T)f{x,T)] 



(5.8) 



which can be used to produce the optimal estimate and 
the associated error of the system process Xr, including 
the magnetic field br = xit- 



B. Smoothing in phase space 

The usual strategy of solving the quantum estimation 
problem is to take the Sx ^ Sy, Sx ^ Sz limit, assume 
Sy and Sz are continuous, and approximate the condi- 
tional quantum state as a Gaussian state [8l.[9L[Tol.[H.[l3|. 
This is akin to approximating the spherical phase space 
with a flat one near S = {Sx,0,0). While the Gaussian 
approximation is probably the most practical, in order 
to illustrate the discrete phase-space formalism proposed 
in Sec. llVi I shall first attempt to convert Eqs. (|5.6p and 
(|5.7p to stochastic equations for discrete Wigner distribu- 
tions in the 2N x 2N phase space before making further 
approximations. 

Let 

m = g-s, N = 2s + 1, (5.9) 

where s is the total spin number. Then 
- _ 2ttp 



N 



(5.10) 



is the operator for the azimuthal angle of the spin vector. 
I shall use m and as the phase-space variables instead 
of q and p, and let 



f{m,4>) ^ f [ q ^ 771 + s,p ^ ^ 

ZTT 



(5.11) 



First consider the measurement-induced decoherence 
term on the second line of Eq. (|5.6p . It can be rewritten 
as 



-iK,q 



/. (5.12) 



Using the definition of the discrete Wigner function given 

by Eqs. (|4.5p and ^ to denote the transform to the 2N x 
2N phase space, it can be shown that 

^m/ ^ ^ J((/) - (l)')f{m, (j)' , x, t) - /(m, (j), X, t), 
<t>' 

(5.13) 



where 



1 r sin[iV(0 -(j)' - k)] 



' m [tan[((/)- 0' - k)/2] 
sin[iV(0 -(/)' + k)] 1 
teLn[{(j) -(/)' + K)/2]j' 



(5.14) 
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While Eq. (|5.13|) has the appearance of the jump term in 
the Chapman-Kolmogorov equation [Eq. (|2.2p ]. J{(f>—(f>'), 
which plays the role of a jump probability density, can 
become negative. In the special case of k = nfi/N, where 
jj, is an integer, however, J{(j> — (j)') is simplified to 



J{(j)~ (j)') = - {54,-n 



(5.15) 



and the measurement-induced decoherence introduces 
random azimuthal jumps in steps of k to the spin vector 
around the z axis. In the limit of N oo, (j) becomes ap- 
proximately continuous, k « 7r/i/iV, and Eq. (|5.13p can 
be rewritten as 



f{m, (j) + K,x,t) + f{m, (j) ~ K,x, t) 
f{m,(j>,x,t). 



(5.16) 



The N oo limit is akin to a ppr oximat ing the spin 
system as a harmonic oscillator [23| and the spherical 
phase space as a cylindrical one. If k <C {A(fP)^/^, we 
can further make the diffusive approximation: 



K 



■J{m,(j),x,t). 



(5.17) 



Next, consider the Larmor 
(ij/h)b[Sy, /]. In terms of rh and 0, 



precession term 



h 

2i 



exp 



(—10) \/s{s + 1) — TO(m + 1) 



exp 



{i(f>)\/ s{s + 1) — rh{rh — 1) 



(5.18) 



With this form, it is difficult to convert the Larmor pre- 
cession term to the phase-space picture analytically, so 
we again make the cylindrical phase-space approxima- 
tion with s ^ (rfi), (Am^)^/^, so that the spin vector 
distribution is concentrated near the equator. This ap- 
proximation is valid when the magnetic field is small and 
fluctuating around zero, or a control, such as an applied 
magnetic field ^, ^, IC^ or an adjustable direction of the 
optical beam, is present to realign the spin vector with 
respect to the optical beam propagation direction. Then 



Sy ~ —hs sin ( 



Cif^ ~jbs cos ( 



(5.19) 

/(m -I- i, (j), X, t) - f{m - i, 0, x, t) 



(5.20) 

Although this looks like the jump term in Eq. (|2.2p . the 
apparent jump probability density is again negative. To 
make the classical connection, assume that m is contin- 
uous and approximate the difference as a derivative: 



^if ^ -jbscos(t)—f{m,(j),x,t), 

which becomes equivalent to the drift term in Eq. 
with Am = jbs cos 4>. 



(5.21) 



Finally, let us consider the terms L±fL!^ in Eq. (|5.6p . 
It is not difficult to show that, in the continuous 6 limit, 



L±fL± ^ |ap <j ^ [/(m, (t) + K,x,t) + f{m, (f> - k, x, i)\ 



± sin(2KTO)/(m, (/), x, t) 



(5.22) 



These terms do not have exact analogs in the correspond- 
ing classical equation [Eq. (|2.1ip ]. unless we make the 
K ^ (AcfP)^/^ approximation. 



L±fLi^\a\^fim,<f>,x,t) + 



■f{m,(l),x,t) 



± sin(2KTO)/(m, (f), x, t) 



\a\^ [1 ± sin(2Km)] /(m, 0, x, t). 



(5.23) 
(5.24) 



Summarizing, a classical model of atomic magnetometry 
can be obtained if we approximate the spherical phase 
space as a cylindrical one near the equator, assume m is 
continuous, and let k <C (Ac^^)-'^/^. The resulting equa- 
tions for /(m, 0, x, t) and g{m, cp, x, t) are 

-t- ^ idN^~dt)iX^-l)f, (5.25) 

-dg = dt ^Chg + lbs cos ^^+^-^^ 

+ ^ idN^-dt)iX^-l)g, (5.26) 



A± = |a|^ [1 ±sin(2Km)] . 
The equivalent system equations are then 



dxt = A{xt,t)dt - 
drrit — dtjbts cos (- 



B{xt,t)dWt., 



(5.27) 

(5.28) 
(5.29) 



where dWt and d4>t are independent Wiener increments 
with dWtdW^ = Q(t)dt and d(j)'^ = \a\'^Kdt. Unlike the 
Gaussian model 0, liol [T^ , this slightly more general 
model shows that the z component of the spin is coupled 
to (j)t via Larmor precession, as one would expect from 
classical dynamics, since Sx ~ hs cos (f> when s ^ m. 
The diffusion of (j> would therefore reduce the estimation 
accuracy in the long run. 

To make the Gaussian approximation, let 
(0t), (A0J )^/^ < 1, so that cos(j)t ~ 1, and let xt 
be a Gaussian random process, such as the Ornstein- 
Uhlenbeck process K K(m), ^(Am^)^/^ < 1, and 

the effective noise covariances are (A±) « |a|^ one can 
use the linear Mayne- Eraser- Potter smoother 0, [3l|, [l^l , 
which combines the estimates and covariances from a 
predictive Kalman filter and a retrodictive Kalman filter, 
to produce the optimal estimate of Xt- Other equivalent 
linear smoothers may also be used (32j . 
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VI. HARDY'S PARADOX IN PHASE SPACE 

In this section, I shall study Hardy's paradox 15] in 
phase space using the quantum smoothing quasiproba- 
bility distribution defined by Eq. (|4.10[) . which allows 
one to estimate quantum degrees of freedom given past 
and future observations in a way closely resembling clas- 
sical estimation theory. The more physical and intu- 
itive Feynman-Wootters distribution is used, since its el- 
ements all correspond to actual paths in the setup. I 
shall show that the paradox arises because the predictive 
Wigner distribution becomes negative, and quantum me- 
chanics becomes incompatible with classical estimation 
as a result. This approach is somewhat different from 
Aharonov et aZ.'s attempt to explain Hardy's paradox 
using weak values [161] . 




FIG. 2: (Color online). Setup of Hardy's paradox. 

As a brief review of the paradox, consider two Mach- 
Zehnder interferometers, one for a positron and one for 
an electron, depicted in Fig. [2l If the interferometers 
are physically separate, then the setup can be configured 
so that the particles always arrive at the C'^ and C~ 
detectors, respectively. Now let's make one arm of an 
interferometer to overlap with an arm of the other. Af- 
ter the first pair of beamsplitters, the two particles may 
meet in the overlapping arms, in which case they annihi- 
late each other with probability 1. With this overlapping 
setup, there is a 1/16 probability that the particles will 
arrive at and , respectively, according to quantum 
theory. 

The paradox arises when one tries to use classical 
reasoning to estimate which arms the particles went 
through. If detects a positron, then the electron must 



have been in the overlapping arm to somehow influence 
the positron to go to D'^ instead of C'^ . The same reason- 
ing can be applied when Z?^ detects an electron, which 
should mean that the positron was in the overlapping 
arm. But if both particles went through the overlapping 
arms, they should have annihilated each other and would 
not have been detected. 

Denote the position of a particle in one arm as and 
that in the other arm as 1 , as shown in Fig. [5J At the 
time instant labeled 0, 



l*)o = |0,0), 



(6.1) 



where the first number in the ket denotes the position of 
the positron, the second number denotes that of the elec- 
tron, and the subscript is the time label. The correspond- 
ing two-particle Wigner distribution using Eqs. (|4.ip and 
is 



/o(0, 0,0,0) /o(0, 0,0,1) 

/o(0, 1,0,0) /o(0, 1,0,1) 

/o(l, 0,0,0) /o(l, 0,0,1) 

/o(l, 1,0,0) /o(l, 1,0,1) 



/o(0,0, 
/o(0,l, 
/o(l,0, 
/o(l,l, 



0) 
0) 
0) 
0) 



/o(0, 
/o(0, 
/o(l, 
/o(l, 



/ 1 1 1 1 



Vo 



0,1,1) 

1,1,1) 

0,1,1) 

1,1,1) 

(6.2) 
(6.3) 



After the first pair of beamsplitters 
1 



l*>i 



(10,0) + 10,1) + 11,0) + 11,1)), 



/ 1 

10 

10 

V 1 



(6.4) 
(6.5) 



If the annihilation did not occur, the a posteriori quan- 
tum state is 



I*) 



1 



V3 



f2{q ,q ,p ,p ) = — 



(10,0) + 10,1) + 11,0)), (6.6) 



(6.7) 




The Wigner distribution has negative elements and can 
no longer be regarded as a classical phase-space proba- 
bility distribution. The negative elements, as one shall 
see later, can be regarded as the culprits that cause the 
paradox. The predictive marginal distributions are still 
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nonnegative, however. In particular, 



//2(0,0) 
/2(0,1) 
/2(1,0) 

V/2(l,l) 




(6.8) 



(6.9) 



which correctly predicts the a posteriori position proba- 
bility distribution if one measures the positions of the 
particles at that instant using strong measurements. 
Most importantly, /2(1, 1) = 0, and the probability that 
one measures both particles in the overlapping arms with 
strong measurements is zero. 

Now perform retrodiction. Given that and 
click, it can be shown that 



52(9^,9 ,P^,P ) 



1 

1 ( 1 

4 I 1 

1 



(6.10) 



The smoothing quasiprobability distribution at time in- 
stant 2 becomes 



Quantum smoothing, on the other hand, is able to 
overrule quantum filtering because some elements of 
72(1, l,p+,_p^) are negative. This way /2(1,1) can 
still be zero with nonzero /2(1, l,p~'',p^) elements, and 
h2{l,l,p~ ,p'^) and /i2(l,l), conditioned upon the de- 
tection outcomes, can become nonzero. The negative 
elements of /2('?^, 9~,P^,P~) thus cause filtering and 
smoothing to produce contradictory trajectories. 

If the detection outcomes are different, say, C+ and 
D~ click, then 



92{q^,q ,P'^,P ) 



h2{q^,q ,P'^,P ) 



/O 1 








10 1 








401 








Vo 1 








/ 

















2 








Vo -1 




















(i) 













(6.14) 



(6.15) 



(6.16) 



h2{q^,q ,p^,p )oc/2((7+,g )g2{q^,q ,p^,p ) 



h2{q^,q ) = 



/O 



Vo 1 

("A 

• 



(6.11) 
(6.12) 

(6.13) 



Hence, given that the annihilation did not occur and 
and D~ click, both particles "reappear" in the overlap- 
ping arms according to quantum smoothing. This result 
is consistent with the classical logic that leads one to the 
same paradoxical conclusion. Mathematically, the para- 
dox arises because the filtering estimation according to 
f2{q^,q~) contradicts the smoothing estimation accord- 
ing to h2{q^ ,q~), with the former ascertaining that the 
particles cannot both be in the overlapping arms, while 
the latter insisting the opposite. 

To see why this cannot happen in classical estimation 



theory, assume for the time being that f2{q^, q ,P^ iP 



) 



is nonnegative. Then /2(1,1) is zero if and only if 
72(1, l,p+,p~) is zero for allp+ andp~. If /2(1, l,p+,p^) 
is zero, so are ft.2(l, 1,P^,P ) and /i2(l,l) according to 
Eq. I|6.1ip . In other words, in classical estimation, if fil- 
tering estimates that the two particles cannot both be 
in the overlapping arms, then no amount of smoothing 
afterwards can alter the certainty of this fact. 



and we have a negative "probability." Leaving aside 
the question of interpreting a negative probability [23 |. 
^2(9^,'? ) still suggests that the most likely positions 
are (g+,(7~) = (1,0), which are consistent with classical 
reasoning and do not contradict the filtering results in- 
dicated by /2 (g"*" , ) . Similarly, when C+ and C~ click, 
the most likely {q^,q^) according to h2{q~,q^) is (0,0), 
which is again what one would expect from a classical 
argument. In this example at least, the most likely po- 
sitions suggested by quantum smoothing coincide with 
the ones obtained by qualitative classical reasoning, as 
depicted in Fig. [31 

To summarize, quantum phase-space filtering and 
smoothing are able to reproduce mathematically the 
salient features of Hardy's paradox and identify the neg- 
ativity of f2{q~^, q~ ,p'^ ,p~) as the culprit that makes the 
classical phase-space picture and quantum theory incom- 
patible. 



VII. CONCLUSION 

In conclusion, the time-symmetric smoothing theory is 
extended to account for discrete variables in classical sys- 
tems, quantum systems, and observations. To illustrate 
the extended theory, atomic magnetometry and Hardy's 
paradox are studied using quantum phase-space smooth- 
ing. The generalized smoothing theory outlined in this 
paper is expected to be useful in future quantum sensing 
and information processing applications. 
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FIG. 3: (color online). The most likely paths undertaken 
by the particles indicated by quantum smoothing given the 
detection outcomes, provided that annihilation did not occur. 
These paths coincide with those suggested by qualitative clas- 
sical reasoning. When and D~ click, the estimated paths, 
as shown in the bottom-right figure, contradict with the fact 
that annihilation did not occur and the two particles could 
not have both been in the overlapping arms. 
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APPENDIX A: OBTAINING THE QUANTUM 
SMOOTHING DISTRIBUTION BY WEAK 
MEASUREMENTS 

In the case of continuous variables, the quantum 
smoothing distribution can be obtained from the statis- 
tics of weak position and momentum measurements, con- 
ditioned upon past and future observations [5] . One may 
also apply a similar method to the discrete- variable case. 
Interestingly, the statistics of weak measurements natu- 
rally lead to a 2N x 2N phase space. 

Consider consecutive q and p measurements of a quan- 
tum system. Let the measurement operators be 



Af-l 



M{yq) = VC ^ exp 

9=0 

M(yp) = x/C ^ exp 



e 27r 

-cos— (y, -g) 



p=0 



e 27r ■ 
2 cos— (yp -p) 



k>(9|, (Al) 
\P){PI (A2) 



where C is a normalization constant and e parameterizes 
the measurement strength and accuracy. The probabil- 
ity distribution of y^ and yp, conditioned upon past and 
future observations, is 

P{yq,yp) 



tr(g/) 



C2 



— } exp 

Ntiigf) , 



e 27r(y, - q) 

- cos 7- 

2 TV 



X exp 
X exp 

X exp 



e 27r(yq - 

- cos \- 

2 N 



e 2iT , , e 27r, ,, 
- cos —(yp - P) + 2 cos j^iVp - P ) 



2Tri{p'q' — pq) 
N 



{p'\g\p){q\fW). (A3) 



Let 

q 



q + q' q' - q _ p + p' 



P -P 



(A4) 



Applying trigonometric identities, one obtains 

27r(y, - q)' 



P{yq,Vp) 



Ntr{gf) 
X exp 



ecos 



J2 

27r(yp -p) 



ecos ■ 



N 



N 



X exp 
X exp 
X exp 



27r(yg ~ q) . 2 ™ 
-2e cos sm — 

N N 

27r(yp -p) . 2 TTU 

-2ecos 7- sm — 

N N 



47ri(wg -|- pu) 



N 



{P + v\g\p -v){q- u\f\q + u). 

(A5) 



Utilizing the periodic nature of the above expression, one 
can change the sum in terms of (g, q') to a sum in terms 
of {q,u), 



N-l N-1 



q=0 q'=a 



q^u 



9e<iO,i,...,iV-i^, 



u e 



N 1 N N 



(A6) 

(A7) 
(A8) 



likewise for {p,p') and {p,v), and the matrix elements 
{p + v\f\p — v) and {q — u\g\q + u) are assumed to be zero 
whenever p + v, p — v, q — u, or q + u is not an integer. 



10 



Thus, 

P{Vq,yp) 

= — > exp 

tr(5/) 

X 9{q,P)f{q,p), 
where 



27r(?/q - q) 27r(yp - p) 

e cos 1- e cos 



TV 



N 



(A9) 



1 

2iV 



^exp 



27r(y5 - g) . 2 t^w 

-2e cos — sm — 

N N 



xexp(^^] {q~u\f\q + u), (AlO) 



~9iq,p) = exp 



27r(yp ~p) . 2 TTW 

-2ecos ^- sm — 

N N 



xexp{^^j{p + v\g\p-v)- (All) 

In the Hmit of infinitcsimally weak measurements and 
e < 1, 

MP) « ^E<^^P (^) (9-"l/l9 + ">, (A12) 



g{q.p) 



which are precisely the discrete Wigner distributions in 
the 2N X 2N phase space. Equation (|A9|) becomes 



P{Vq.Vp) 



exp 



9,P 



ecos 



iV 



TV 



Hq,p), 

(A14) 



and can be regarded, from the perspective of classical 
probabihty theory, as the probabihty distribution for 
noisy q and p measurements, when the system has a 
phase-space distribution given by the quantum smooth- 
ing distribution h(q,p). h{q,p) can therefore be obtained 
in an experiment with smaU e by measuring P{yq, Vp) for 
the same g and / and deconvolving Eq. (|A14p . 
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